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ABSTRACT 

There is mounting observational evidence from Chandra for strong interaction between keV gas and 
AGN in cooling flows. It is now widely accepted that the temperatures of cluster cores are maintained 
at a level of ~ 1 keV and that the mass deposition rates are lower than earlier RO SAT /Einstein values. 
Recent theoretical results suggest that thermal conduction can be very efficient even in magnetized 
plasmas. Motivated by these discoveries, we consider a "double heating model" which incorporates the 
effects of simultaneous heating by both the central AGN and thermal conduction from the hot outer layers 
of clusters. Using hydrodynamical simulations, we demonstrate that there exists a family of solutions 
that does not suffer from the cooling catastrophe. In these cases, clusters relax to a stable final state, 
which is characterized by minimum temperatures of order 1 keV and density and temperature profiles 
consistent with observations. Moreover, the accretion rates are much reduced, thereby reducing the need 
for excessive mass deposition rates required by the standard cooling flow models. 

Subject headings: galaxies: clusters — cooling flows — X-rays: galaxies — conduction — intergalactic 
medium 



1. INTRODUCTION 

Radiative cooling of gas in the central regions of galaxy 
clusters occurs on a timescale much shorter than the Hub- 
ble time. The cooling time-scale of this gas increases 
with distance from the cluster core. In the absence of 
any heating sources, this implies that the intracluster 
medium must accrete subsonically toward the center in 
order to maintain pressure equilibrium with the gas at 
larger radii. The mass deposition rates predicted by 
this cooling flow model are very high and range typi- 
cally from 10 to 1000 solar masses per year. X-ray ob- 
servations made prior to the launch of the Chandra Ob- 
servatory seemed to be broadly consistent with this pic- 
ture (Fabian 1994). Although both the gas temperature 
and cooling time are observed to decline toward cluster 
cores, new Chandra and XMM-Newton observations show 
a remarkable lack of emission lines from gas at temper- 
atures below ~ 1 keV in the central regions of clusters 
(Peterson et al. 2001; Allen et al. 2001). Moreover, the 
cooling mass deposition rates obtained with Chandra and 
XMM using spectroscopic methods are ^ 10 times smaller 
than earlier estimates based on ROSAT and Einstein 
observations (McNamara et al. 2001; David et al. 2001; 
Peterson et al. 2001). On the other hand, morphological 
cooling rates give accretion rates vastly exceeding the ones 
based on spectroscopic methods (David et al. 2001). The 
strong discrepancy between these results indicates either 
that the gas is prevented from cooling by some heating 
process or that it cools without any spectroscopic signa- 
tures (Fabian 2001). 

In this paper, we consider cooling flow models that in- 
corporate the effects of heating by central active galactic 
nuclei and thermal conduction from the hot outer layers of 
clusters. We show that evolving density and temperature 
profiles can relax to stable final states, which are consis- 
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tent with X-ray observations. Equilibria are characterized 
by minimum temperatures of order 1 keV at radii ~ 1 kpc. 
The paper is organized as follows. In section 2 we discuss 
the observational and theoretical motivation for including 
both AGN heating and thermal conduction. In section 3 
we present the details of our model. The results of hydro- 
dynamic simulations are presented in Section 4 and the 
main results are summarized in section 5. 

2. HEATING MECHANISMS 

2.1. Feedback and heating by ACN 

A large fraction (~ 70%) of cD galaxies in the centers 
of cooling flow clusters shows evidence for powerful radio 
activity. This fraction is much larger than in non-cooling 
flow clusters and suggests a link between the presence of 
the cooling gas and the activity of the central supermas- 
sive black hole. Currently, Chandra offers the best op- 
portunity to study this coupling as the cluster cores are 
optically thin and can be spatially resolved down to scales 
of a few kpc within z < 0.1. These scales are comparable 
to the sizes of the radio sources. Recent Chandra obser- 
vations of the Perseus, Hydra A and many other clusters 
reveal holes in the X-ray surface brightness which coincide 
with radio lobes (Fabian 2001; McNamara et al. 2001). 
These radio sources inflate bubbles of hot plasma that 
subsequently rise through the cluster atmosphere, mix- 
ing and heating it up. This mechanism has been invoked 
to explain the reduced cooling rates (Jones et al. 2002; 
Churazov et al. 2002). 

2.2. Heating by thermal conduction 

A potential difficulty with reducing the cooling accre- 
tion rate by means of a central heating source is that too 
strong a heating source will lead to very strong convec- 
tion, which could remove metallicity gradients observed 
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by Chandra in cluster cores (Johnstone at al. 2002). 
Therefore, the hmits on energy requirements of cen- 
tral AGN and on the amount of mixing and convec- 
tion suggest that an additional heating source is re- 
quired to offset cooling rates. An important candi- 
date is thermal conduction of heat from the outer hot 
layers of the cooling flow cluster. With some excep- 
tions (Bertschinger and Meiksin 1986; Malyshkin 2001; 
Brighenti and Mathews 2002; Voigt et al. 2002; 

Fabian ct al. 2002) this effect has been largely neglected 
in observational and theoretical studies due to the pre- 
conception that the conduction coefficient would be sup- 
pressed much below the classical Spitzer value by mag- 
netic fields. Indeed, there is strong evidence that clus- 
ters are magnetized by intermittent radio galaxies in 
cluster cores (Clarke et al. 2001; Kronberg et al. 2001; 
Reynolds and Begelman 1997). However, recent theoreti- 
cal work by (Narayan and Mcdvedcv 2001) suggests that 
if the magnetic field is highly turbulent, then thermal 
conduction is relatively eflacient. They find that if the 
turbulence spectrum extends over two or more decades 
in wave vector, then the thermal conduction coeSicient is 
only a factor ^ a few below the Spitzer value and thus 
could play a significant role in cooling flows in clusters 
of galaxies. The other argument in favor of thermal con- 
duction is that it is a very strongly increasing function of 
temperature and, therefore, could be efficient in the outer 
layers of clusters where temperatures are high. 



2.3. Simultaneous heating by AGN and conduction 

If radiative cooling is balanced solely by energy input 
from a central supermassivc black hole then a cooling 
flow can be quenched provided that a sufficient amount of 
power is provided in kinetic form (McNamara 2002). Mod- 
els based on the assumption that radiative cooling is bal- 
anced by energy input from AGN predict that the mechan- 
ical power of AGN in cooling flows is much higher than the 
presently observed bolometric luminosities of these objects 
(Churazov et al. 2002). However, the power necessary to 
offset a cooling flow may still be present in the kinetic 
form after the central AGN becomes inactive. In fact, the 
required power is consistent with estimates of jet power 
in some objects. Theoretical models of impulsive central 
heating in elliptical galaxies (Binncy and Tabor 1995) pre- 
dict violent successive cooling catastrophes and tempera- 
ture rising towards the center for radii less than ~ 100 kpc. 
On the other hand, attempts to build cooling flow models 
with conduction as the only heating mechanism fail. Such 
models predict that the cooling time of the central clus- 
ter region diminishes to less than the free-fall time, lead- 
ing to a cooling catastrophe (Meiksin 1988). This results 
in supersonic accretion of cold and dense material in the 
core, which is contrary to observations. Stable models can 
be obtained but they assume a significant amount of dis- 
tributed mass dropout which does not have observational 
grounding. Moreover, even when the models are stable, 
such cooling fiows suffer from additional problems. Due 
to the strong dependence of the conduction coefficient on 
temperature (cx T^^^), the predicted temperatures either 
drop by less than a factor of 2 or the central tempera- 
tures are so low that heat conduction is negligible and 



the central accretion rate often becomes very high. This 
suggests that heat conduction alone cannot substantially 
reduce mass accretion rates everywhere within a cluster. 
Recently, (Brighenti and Mathews 2002) considered cen- 
trally heated cooling flow models with conduction. They 
concluded that such models are grossly incompatible with 
observations and were unable to find even a single accept- 
able cooling flow model using their heating prescriptions. 

In the next section we present details of our model and 
look at the consequences of bridging these two heating 
regimes. 

3. "double heating" model 

3.1. Physical assumptions 

We solve numerically the equations of hydrodynamics in 
the following form 



dp 

di 

dS 

'dt 

de 

dt 



+ V • (pv) = 
+ V • (Sv) = -Vp - pV* 

V • (ev) = -pV • v - V • Fcond 

-nlA{T) + H, 



V • F 

V J- n 
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where p is the mass density, rie the electron number den- 
sity, p the pressure, v the velocity, S = pv the momen- 
tum vector, e the internal energy density, "if the gravi- 
tational potential, A(T) the cooling function and H the 
heating rate per unit volume. We adopt an equation of 
state p = {■J — l)e and consider models with 7 = 5/3. The 
conductive flux Fcond is given by 



Fcond = -/kT5/2VT, 

where k is the Spitzer conductivity 



K = 



1.84 X io-'"^r5/2 

In A ' 



(4) 



(5) 



with the Coulomb logarithm In A = 37, and / is a reduc- 
tion factor (0 < / < 1). The saturation of the conductive 
flux was not important for the parameters considered in 
this paper. The convective fiux Fconv is given by mixing 
length theory 



conv — 1 Q 



V/V™(-Vs)3/2 ifV.i<0 



otherwise 



(6) 



where g is the gravitational acceleration, is the mix- 
ing length, s ~ Cy lii[(7 — l)e//9^] is the gas entropy, Cy is 
the specific heat per unit volume, Cp ~ ^kBl\(ci ^ 1)^*^^] 
is the specific heat per unit mass at constant pressure, 
p, = 0.5 is the mean molecular weight and 7 is the adi- 
abatic index. In mixing length theory Im is a free pa- 
rameter. We use Z„i ~ min[0.3(7 — l)e/{pg),r], where 
r is the distance from the cluster center. Motivated 
by recent Chandra X-ray observations of clusters (e.g., 
(Schmidt et al. 2001; Allen et al. 2001)), we parameterize 
the dark matter distribution using a non-singular isother- 
mal mass density profile. We add a contribution from the 
central galaxy also in the form of a non-singular isothermal 
profile but characterized by different values of parameters. 
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We assume that the total potential does not evolve. The 
contribution of the gas to gravity is negligible throughout 
the simulation volume. The corresponding total gravita- 
tional potential is then given by \E'tot = *I'c + ^g, where 



1 + 



(^) -ctan 



(7) 

where ac,g is the cluster or galactic velocity dispersion and 
rc,g is the core radius of the cluster or galaxy. 

Note that equations (l)-(3) do not include terms related 
to mass dropout rate. Such terms arc usually invoked 
to avoid a cooling catastrophe by removing cooled gas 
from the flow. Thermal instabilities, which lead to mass 
dropout, appear when the cooling time is comparable to 
or shorter than the dynamical time. In our simulations the 
cooling time is longer by a factor 10—100 than the dynam- 
ical time, so our assumption is justified. A small amount 
of mass dropout may still occur in the central regions if 
the gas becomes locally ovcrdcnsc, e.g., due to interactions 
of the central AGN with the surrounding medium. How- 
ever, as mentioned above, there is no strong observational 
indication that significant amounts of gas are decoupling 
throughout the cooling flow region, although star-forming 
regions and small amounts of cool gas in the form of Ha 
emission have been detected in the innermost regions of 
clusters. Such colder gas also could have been lifted from 
very small radii to larger distances by the central AGN. 



3.1.1. Heating and cooling 

Following (Tozzi and Norman 2001) we use an approxi- 
mation to the cooling function based on detailed calcula- 
tions by (Sutherland and Dopita 1993) 



nlK = [CiikeTY + C2{kBTf + Cs]n,n, 



(8) 



where is the ion number density and the units for 
ksT are keV. For an average metallicity Z = 0.3^0 
the constants in equation (8) are a = —1.7, /3 = 0.5, 
Ci = 8.6 X 10^^ C2 = 5.8 X IQ-^ and C3 = 6.4 x 10"^ and 
we can approximate riirie = 0.704(p/mp)^. The units of A 
are 10~^^ erg cm"^ s~^. This cooling function incorporates 
the effects of free-free and line cooling. 

The final term in equation (3) represents heating by the 
central AGN. Since it is ultimately heating by mechani- 
cal energy and particles from radio sources, it will be dis- 
tributed in radius. We adopt a physically motivated "effer- 
vescent heating" mechanism (Begelman 2001) to describe 
the volume heating rate. The physical motivation for this 
model is as follows. Suppose the central radio source de- 
posits buoyant gas which distributes itself relatively evenly 
among bubbles or filaments but does not mix microscopi- 
cally with the intracluster medium (ICM). These bubbles 
will then rise through the ICM and expand because of the 
non-negligible pressure gradient. Bubble expansion will 
be associated with the conversion of the internal bubble 
energy to kinetic form and, eventually, to heat due to dis- 
organized motion of the ICM. Since the bubbles rise on 
a timescale shorter than the radiative cooling time, the 
mechanism should reach a quasi-steady state and the de- 
pendence on details of bubble filling factor, rise rate, etc. 



cancel out. In a steady state (and assuming spherical sym- 
metry), the energy flux available for heating is then 

eocp6(r)('^''-i)/T\ (9) 

where Pbir) is the partial pressure of buoyant gas inside 
bubbles at radius r and 75 is the adiabatic index of buoy- 
ant gas. Assuming that the partial pressure scales with the 
thermal pressure of the ICM, the volume heating function 
Ti can be expressed as 



n 



-/i(r)V- 



(76-l)/7!> 



1 dlnp 



, (10) 



Anr'^ ' \po J r dlnr' 

where po is the central pressure and h{r) is the normaliza- 
tion function 



h{r) 



(l_e-'-/'-o)g-i. 



(11) 



In equation (11), L is the luminosity of the central source 
and q is defined by 



q 



r 



PoJ 



(7b-l)/7t 



1 dlnp 



(1 



,-r/r, 



)dr, (12) 



r dlnr 

where ro is the inner heating cutoff radius. Interest- 
ingly, the cosmic ray heating model for cooling flows 
(Loewenstein et al. 1991) predicts a very similar func- 
tional form for the heating function. In this model, heating 
is due to cosmic rays, which are produced in the central ac- 
tive nucleus and are trapped by Alfven waves in the cooling 
flow plasma. They also do work on the thermal gas as they 
propagate down the pressure gradient. Radio galaxies are 
likely to be intermittent on a time scale much shorter than 
the Hubble time, and possibly as short as ti ^ 10^ — 10'"' 
yr (Reynolds and Begelman 1997). The bubble rise speed 
would be comparable to the sound speed, which in turn is 
comparable to the dynamical speed. Therefore, the cool- 
ing timescale is much longer than the bubble rise timescale, 
and it is justifiable to treat feedback heating i in a time- 
averaged sense. Therefore, in equation (11) we assume 
that the luminosity is injected instantaneously and ne- 
glect any delay between central activity and heating of 
the ICM. We also assume that all energy generated in the 
center goes into heating. In principle some fraction of the 
energy could escape the cluster in the form of sound waves. 
However, they are likely to carry away only a fraction of 
the energy comparable to that which is dissipated (e.g., 
(Reynolds et al. 2002; Churazov et al. 2002)). The phys- 
ical motivation for the inner heating cutoff arises from the 
finite size of the central radio source. Energy from such 
a source will start to dissipate at a finite distance from 
the cluster center. Unlike heating functions that depend 
on local microscopic physics, this heating function is non- 
local in the sense that it depends on the pressure gradi- 
ent. In this regard, it resembles thermal conduction, but 
here the heating rate depends on the gradient of pressure 
rather than temperature. For simplicity and in order to 
avoid numerical complications, we adopt hydrostatic val- 
ues to calculate the logarithmic derivative in the code (i.e., 
Vp= -pV*). 

In the specific case considered below, the feedback lumi- 
nosity is assumed to be L ~ —eM(?, where M = Anr^^^^pv 
is the accretion rate at the inner radius rmin and e is the 
accretion efliciency. 
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3.2. Numerical methods 

All of the calculations presented in this paper use 
the ZEUS-3D (Clarke ct al. 1994) code in its ID mode. 
The code has been modified by including a non-singular 
isothermal potential, feedback heating, cooling, convection 
and conduction. For stability, the convection and conduc- 
tion terms have to be integrated using time steps that 
satisfy appropriate Courant conditions. The Courant con- 
ditions for conduction and convection read 



At, 



^ 1 e(Ar) 



2 /kT5/2 



125^(Ar)^ 



(13) 

(14) 



Whenever either (or both) of these timescales is smaller 
than the time step used for hydrodynamic equations, we 
integrate the respective terms separately at the smaller 
time step. This method is analogous to the one used by 
(Stone et al. 1999) in hydrodynamical simulations of vis- 
cous non-radiative accretion discs. 

Our computational grid extends from rmin = 1 kpc to 
^■max = 200 kpc. In order to resolve adequately the inner 
regions it is necessary to adopt a non-uniform grid. We 
use a logarithmic grid in which (Ar)i+i/(Ar)i = "~v^, 
where A'' is the number of grid points. Our standard res- 
olution is N = 400. 

3.2.1. Initial and boundary conditions 

In order to generate initial conditions we solve the equa- 
tion of hydrostatic equilibrium of a gas at constant temper- 
ature. We assume that the gas density initially constitutes 
a fraction /igas of the cluster dark matter density at the 
inner radius. The gas is assumed to be in contact with a 
pressure and thermal bath at the outer radius. Thus, we 
ensure that temperature and density at the outer radius 
are constant. We adopt inflow/outflow boundary condi- 
tions at the outer radius but at the inner boundary wc use 
outflow boundary conditions. We extrapolate hydrody- 
namic variables from the active zones to the "ghost" zones 
used to compute the derivatives of the hydrodynamic vari- 
ables. At the inner boundary we allow mass to flow out 
of the computational region {v < 0) and use a switch to 
ensure that the fluid can only flow toward the center (i.e., 
boundary value of the velocity is zero if the gas tries to 
flow in the opposite direction). Feedback is assumed to 
be present only when the gas at the inner boundary flows 
inwards {v <0). 

4. RESULTS 

To illustrate our method we now present one repre- 
sentative model. This model reproduces the main fea- 
tures of the observed cooling flows including a floor in 
the temperature at about 1 keV. We note that fine tun- 
ing is not required to obtain stable equilibrium solutions. 
Figure 1 shows the evolution of temperature (upper left 
panel), electron number density (upper right panel), en- 
tropy (S = ksT /n1^^] lower left panel) and mass ac- 
cretion rate as a function of distance from the cluster 
center. We follow cluster evolution for one Hubble time 
tH = {Ho = 75 kms-^Mpc"^). Each curve within a 



given panel is plotted every 1/40 or 1/100 of the Hubble 
time, as specified in the caption. The initial gas tempera- 
ture is set to be 4 keV throughout the cluster. Specific pa- 
rameters of the model presented in Figure 1 are as follows: 
e = 0.003, 7fc = 4/3, 7 = 5/3, ro = 20 kpc, r^ = 30 kpc, 
rg = 4 kpc, CTg = 800 km s'^, = 300 km s'^, f = 0.23, 



0.3Zq and figas 
the final gas mass fraction /gas('' < 200 kpc) = Mgas(< 



0.025 (which corresponds to 

r)/Mtot(< r) ^ 0.06). 

Line and free-free cooling in the cluster center leads to 
a slow decrease in temperature. At this stage feedback is 
not yet very important. The initial phase of slow cooling 
is followed by a gradual increase in cooling rate, which is 
caused mainly by the increase of density in the center. The 
feedback heating is controlled by the accretion rate in the 
center. Therefore, the cluster does not cool in runaway 
fashion. Once the average inflow velocity in the center has 
become sufficiently high and sufficient density has accu- 
mulated in the center of the cluster, the feedback becomes 
strong enough to suppress further decrease in tempera- 
ture. The cluster then quickly relaxes to an equilibrium 
state. The evolution of gas density is similar. Initially, 
gas accumulates gradually in the center of the cluster. 
When feedback becomes strong, further accumulation is 
prevented and the central density is stabilized. Entropy 
and mass accretion rate exhibit a similar behavior - the 
former relaxes to a stable equilibrium and the latter, after 
oscillations through positive and negative values, tends to 
a small constant negative value as a function of radius. 
In the final state, the accretion rate at the inner radius 
is M ^ — 1.76Af0 yr~^, which is much smaller than acrc- 
tion rates inferred from the standard cooling flow models. 
Our M corresponds to the bolometric feedback luminos- 
ity L ^ 3 x 10''^ erg s""'^. The most important aspect of 
the chistcr evolution is that it is possible for the cluster 
to reach sustainable and stable equilibrium. In our model, 
there is no need for mass dropout distributed through- 
out the cooling flow region. However, some fraction of 
the material in the central regions may become thermally 
unstable and form Ha emitting filaments or stars (e.g., 
(Blanton et al. 2001; McNamara 2002)). Adding mass 
dropout terms corresponding to these effects would only 
improve the stability properties of our models. It could 
also decrease the central gas density and, thus, lead to a 
higher minimum temperature. 

Figure 2 shows the temperature, electron number den- 
sity, entropy and mass accretion rate as a function of 
time (in units of Hubble time tn = Hq^, Hq = 75 km 
s^^ Mpc""'^) for different distances from the cluster cen- 
ter r = 5, 10, 20, 40, 80, 160 kpc. The parameters 
are the same as for Figure 1, i.e.. Figure 2 shows cross- 
sections through the panels in Figure 1 at the above radii. 
In the case of entropy and temperature, the above se- 
quence of r corresponds to the curves from bottom to top. 
For density the trend is opposite. In the case of accre- 
tion rate, the amplitude of oscillations increases with r 
(= 5, 10, 20, 40, 80). The oscillations are due to sound 
waves, which propagate across the cluster as it adjusts 
to changing conditions. The precise character of these 
sound waves, i.e., number of peaks per unit time, depends 
on the boundary conditions. As can be clearly seen, the 
slow initial evolution of the cluster is followed by a faster 
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"collapse" stage around t = 0.2tH- After this phase, the 
cluster stabilizes and reaches its final equilibrium state. 
The profiles shown in Figure f bear close resemblance 
to the observed temperature, density and entropy pro- 
files (e.g., (Johnstone at al. 2002) Chandra observations 
of Abell 2199). 

Although nominally present in the code, convection is 
not important for the parameters of the model presented 
in this paper. The intracluster medium usually does not 
become convcctivcly unstable. We also considered mod- 
els with stronger heating in the center (i.e., higher effi- 
ciency e). In such models the gradient of entropy was 
significantly negative and convection was present. How- 
ever, strong heating prevented the gas from accreting in 
the center. This led to the accumulation of material and, 
subsequently, to a cooling catastrophe. 

5. CONCLUSIONS 

We have proposed a new class of time-dependent cooling 
flow models where cooling is offset by a combination of cen- 



tral AGN heating and thermal conduction from the outer 
regions. Our models do not require any mass dropout rate 
distributed throughout the cluster. We showed that it is 
possible to obtain stable flnal equilibrium states, which 
do not suffer from the cooling catastrophe. We have pre- 
sented a representative model, which reproduces the main 
features of observed cooling flows, including a floor in the 
temperature at about 1 keV. We have also found stable 
models for other parameters, which we will present later. 
Fine tuning is apparently not required to obtain stable 
equilibrium solutions. Moreover, stable models are char- 
acterized by gas accretion rates that are much smaller 
than the mass dropout rates predicted by standard cool- 
ing flow models. 

We are grateful to Phil Armitage, Fabian Heitsch, Chris- 
tian Kaiser and Daniel Proga for helpful discussions and 

the referee for a fast response. This work was supported 
in part by NSF grant AST-9876887. 
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Fig. 1. — Time sequence of temperature (upper left panel), electron number density (upper right), entropy (5 = fesT/ng ^ , lower left) and 
accretion rate. Temperature and density are shown every O.Olif^'^ and entropy and accretion rate every 0.025^^^^^. The model settles down 
to a stable equilibrium state, which is visible via the dense concentration of curves. See text for additional information. 




Fig. 2. — The dependence of temperature (upper left panel), electron number density (upper right), entropy (lower left) and accretion rate 
as a function of time for different distances from the cluster center. The final accretion rate is ~ 1.76Mq yr~^. See text for details. 



